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Abstract. Two recent numerical studies of hypersonic flows over bodies with 
concavities revealed problems with convergence to a steady state with an oft-used 
application of local-time-stepping. Both simulated flows showed a time-like, periodic 
shedding of vortices in a subsonic domain bounded by supersonic external flow 
although the simulations, using local-time-stepping, were not time accurate. Simple 
modifications to the numerical algorithm were implemented to enable implicit, first- 
order accurate in time simulations. Subsequent time-accurate simulations of the two 
test problems converged to a steady state. The baseline algorithm and modifications for 
temporal accuracy are described. The requirement for sub-iterations to achieve 
convergence is demonstrated. Failure to achieve convergence without time accuracy is 
conjectured to arise from temporal errors being continuously refocused into a subsonic 
domain. 
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1 INTRODUCTION 

Two recent numerical studies of hypersonic flows over bodies with concavities 
revealed problems with convergence to a steady state with an oft-used application of 
local-time-stepping. A concavity is defined as any region over which outward running 
normals to the surface intersect. Compression corners are included in this definition of 
concavity. 

The first study 1 involved a code validation experiment 2 of Mach 9.5 flow over a 
shaip, double-cone at Re/m = 264000 (Run 24) in which experimental data indicated 
steady flow but the numerical simulation failed to converge. In contrast, a numerical 
simulation of Run 28 over the same configuration at Mach 9.6 but with a lower 
Reynolds number (Re/m =144000) converged to a steady state for all grids and physical 
models tested. At the time, it was thought that the non-linear Symmetric Total Variation 
( STVD ) 3 limiter caused the problem because an unlimited, first-order simulation 
converged (albeit a poor simulation of the experiment.) 

The second study involved towed, inflatable, ballute (balloon-parachute) decelerators 
in which hypersonic flow over concave shapes would not support a steady, external 
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hypersonic flow 4 . A spherical segment depression with rounded shoulders supported a 
pulsing wave in which the bow shock oscillated over a distance exceeding the outer 
radius of the ballute. (The simulation was for a uniform, incoming flow - it did not 
include the wake of the towing spacecraft.) The existence of an unsteady flow in this 
case was not suiprising; cavities in the stagnation region of blunt bodies are known to 
induce unsteadiness 5 . 

Both simulated flows showed a time-like, periodic shedding of vortices in a subsonic 
domain bounded by supersonic external flow although the simulation, using local-time- 
stepping, was not time accurate. Simple modifications to the numerical algorithm were 
implemented to enable implicit, first-order in time simulations. Subsequent time- 
accurate simulations of the two test problems converged to a steady state. The baseline 
algorithm and modifications for temporal accuracy are described. The requirement for 
sub-iterations to achieve convergence is demonstrated. Failure to achieve convergence 
without time accuracy is conjectured to arise from temporal errors being continuously 
refocused into a subsonic domain. The failure to converge using local-time-stepping is 
not evident in all simulations of hypersonic flows over concavities; the dependence on 
specific details of the algorithm, geometry, and flow parameters require additional 
study. 

All of the numerical simulations discussed herein are executed using the Langley 
Aerothermodynamic Upwind Relaxation Algorithm (LAURA). 6,7 This paper reviews 
implicit algorithms as applied in LAURA based on single-level storage local-time- 
stepping and two-level storage using sub-iterations to achieve first-order time accuracy. 
These details are reviewed because they strongly influence the evolution of an 
“unsteady” response in the test cases. 

2 ALGORITHM 

2.1 Point-Implicit Relaxation 

The Langley Aerothermodynamic Upwind Relaxation Algorithm (LAURA) employs 
two relaxation modes: point-implicit and line-implicit. The point- implicit mode 
implicitly_solves only contributions of dependent variables at cell L, q L , to the residual 
at cell L, r L ; contributions from dependent variables at neighbor cells are treated 
explicitly. The contributions of inviscid, viscous, and source terms are included in the 
point-implicit matrix, Bl, as follows: 

di L = C 2*^L ( 1 ) 

where Sl is the Jacobian of the chemical and thermal source term components of the 
residual, Ml is the Jacobian_of the inviscid terms, and N L is the Jacobian of the viscous 
terms - all with respect to q L . Scalar constants ci> 1 .5 and C2>0.5 are relaxation factors. 
Typical values are ci=3 and C 2 =l. Relaxation factors on the source term Jacobian have 
never been required. 

Within the context of flux difference splitting (FDS) and Roe’s averaging 8 the 
inviscid Jacobian is defined: 


A = T X I A m | (2) 

^ m 

where m is summed over all faces of cell L, A m is the Roe’s averaged flux Jacobian at 
face m, and |A m | is defined using the absolute values of the eigenvalues of A m . The 
elements of A m are required in the definition of the FDS formulation of the inviscid 
flux; consequently, it is economical to reuse these quantities in the formulation of the 
Jacobian. The absolute value of the eigenvalues may be limited from below using a 
variation of Flarten’s entropy fix 9 "’. 


2 
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The viscous Jacobian is defined: 

A = X 

m 


( 3 ) 


where B m is the Jacobian of viscous terms across the m th face of cell L. It implicitly 
treats only the divided difference across the cell face. Transport properties at cell face m 
are linear averages of adjacent cell center values. Velocities in the shear work term are 
the Roe averaged quantities already computed for the inviscid flux function. 

Eigenvalues of B m are positive or zero. 

The change in dependent variables after a single, point-implicit relaxation step is 
computed as follows: 


-±- I + 

|>t L 

where I - the identity matrix, Q L - the volume, and At L - the time step at cell L (usually 
derived from a constant Courant number, local-time-stepping, specification) are 
associated with the time derivative term of the conservation equations. In the present 
study, At L is based on Courant numbers of 10 1 and 10 6 as well as constant values for all 
indices. 


a 4l = V 


(4) 


2.2 Line-Implicit Relaxation 

The point-implicit option is not efficient for relaxing the viscous layers. Some 
efficiency can be recovered by using multitasking with majority of cycles in the 
boundary layer or by using an over-relaxation factor 1>C2>0.5. However, a better 
approach is to introduce implicit dependence of the solution across the entire shock 
layer. 

The line-implicit relaxation algorithm is a simple extension of the point-implicit 
algorithm. However, approximations to A L , the Jacobian of r L with respect to q L _j , 

and Cl, the Jacobian of ? L with respect to q L+1 , are introduced to preserve symmetry of 
form with Bl that also improves convergence. They are defined: 


<k = 

-f 

+ A L-1 / 2 ) C 2^L-1 / 2 

^L-l / 2 

(5a) 

^ = 


| _ ^L+l/ 2 ) _ C 2 B L+l/2 

^L + l / 2 

(5b) 


It is assumed that sequential indices L define a continuous series of cells spanning some 
portion of the computational domain. The subscripts L ± 1 / 2 refer to Roe averaged 
conditions on face m between cell centers L and L ± 1 , respectively. A more exact 
linearization of the residual would lead to use of Al+/- i in place of Al+/. 1/2 in the 
second occurrence of Eqs. 5; however, the present formulation provides exact 
cancellation of convective, downwind implicit influences and exhibits more robust 
convergence characteristics, particularly in the vicinity of strong shock waves at large 
Courant numbers. 

In the present work, the source term requires no off-diagonal, implicit conmbution. 
The line-implicit relaxation in Eq. (6) requires approximately 2.8 times more CPU time 
per iteration as the point-implicit algorithm in the applications described here but 
convergence rate more than makes up for this overhead. It requires additional storage 
for the off-diagonal blocks. 


3 
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2.3 Implicit Boundary Conditions 

Implicit treatment of boundary conditions is required to obtain good convergence 
with line relaxation. Boundary conditions in the present formulation are implemented 
with ghost cells. If the off-diagonal Jacobian at the ghost cell for L=1 is Ao then Bi is 
augmented by Ao Wo' l W\ where: 


W 0 dq 0 = ‘Mjdq, (7) 

is an expression of the dependence of the ghost cell on the boundary cell. In like 
manner, Bn is augmented by Cn+iI^n+i ^n where W x+l dq N+l = W s dq N . 


2.4 Time Accuracy 

The LAURA algorithm, in its default mode, uses single level storage and relaxes the 
steady form of the governing equations with a Courant number of 1 0 6 . However, some 
applications considered here require time-accurate relaxation. Time accuracy (first- 
order in At) can be recovered by using constant time step advancement at all cells and 
by employing two levels of storage with sub-iterations to converge the current iterate 
before advancing to the next time level 10 . 


At 


I + 



-k.k + 1 
r L 


( 8 ) 


The modified algorithm is conceptually represented by Eq. (8), identical to Eq. (4) 
except that sub-iterations on Aq k are accommodated and a constant time step is used. 
Equation (6) is modified in a similar manner. The superscript (k, k+1) signifies that the 
latest available iterate on q n + Aq k is utilized in the evaluation of r L . Also, r L now 
includes the time derivative term (q} , + l - q“)Q L / At • In theory, Eq. (8) is iterated on all 
cells until Aq k+1 = q k+I - q k is converged to order At 2 for first-order temporal 
accuracy. In practice q n+I is updated after 10 to 20 sub-iterations on index k before 
advancing to time step n+2. 

3 APPLICATIONS 

Two configurations are discussed in which a steady state (or near steady state) could 
only be computed by executing a time accurate simulation. The steady state solution is 
believed to be correct for the ideal boundary conditions applied here in which all 
boundary conditions are steady and the incoming flow is uniform. Possibly, these 
geometries are also more sensitive to flow non-uniformities that trigger unsteady 
response in ground-based or flight experiments. 


4 
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3.1 Concave Ballute 

1 1 12 

Large, inflatable ballutes (balloon-parachutes ) have been proposed ’ as hypersonic 
decelerators for planetary aerocapture applications. Simulations of candidate 
configurations (spheres, toroids) for towed decelerators have been presented previously. 
Parachute-like configurations were also simulated in the early phases of the study but 
were abandoned because of concerns regarding flow stability and aerodynamic stability. 
The chute configurations failed to converge using the constant Courant number 
relaxation scheme. Rather, they showed a periodic growth and collapse of the shock 
wave, the extremes of which are shown in Figure 1 for a simple Mach 6, perfect gas 
flow. Instabilities have been reported previously in blunt body flows with cavities 5 . 
While the magnitude of these chute instabilities was larger than previously reported for 
the blunt body with cavity, their occurrence did not seem unreasonable at the time. 

Additional variations of the relaxation scheme have been tested on this case - all 
starting from the same initial condition. The observations that follow were taken from 
movies in which a new frame is shown every 10 th time step. The solution was relaxed 
using the baseline, local-time stepping, point-implicit algorithm with Courant numbers 
of 10°, 10 1 , and 10 6 . All cases showed unsteady, periodic motion. The amplitude of the 
shock motion was largest for the largest Courant number. 

This test was repeated using the constant time step but varying the number of sub- 
iterations over values of 1, 10, and 20. Larger number of sub-iterations promotes 
improved time accuracy. Point-implicit relaxation is used unless otherwise noted. The 
unsteady response is still very strong with only one sub-iteration (baseline algorithm 
except with constant time step rather than constant Courant number.) The test with ten 
sub-iterations exhibited small but persistent oscillations that never damped over 25 
cycles. The periodic response with constant amplitude was observed after 12 cycles. 

The test case with 20 sub-iterations achieved a near steady state ( shock stand-off varied 
less than 1% after approximately 8 cycles. ) Repeating this case with line-relaxation 
from the same initial condition led to a steady state in fewer relaxation steps. Line 
relaxation provides a greater level of convergence within the sub-iterations - 
consequently the temporal error, while still only first-order accurate, is reduced. 

The effect of sub-iterations on the convergence of a constant Courant number 
simulation was studied. Twenty sub-iterations were executed on a case in which the 
Courant number was fixed at 10 1 . This simulation also failed to achieve a steady state. 



Figure 1 : Density contours over chute at Mach 6 showing unsteady range of motion using constant 
Courant number simulation and steady solution from time -accurate simulation. 
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The LAURA code defaults to use of an alternating directional sweep through the 
domain from the body to the inflow boundary and back. Changing the sweep direction 
( axis to outflow and back) did nothing in the constant Courant number tests to promote 
convergence to a steady state. 

In summary, the test problem requires a time accurate simulation to evolve to a 
steady state. Relaxation using a constant time step but insufficient sub-iterations to 
attain first-order time accuracy fails to produce a steady state. Multiple sub-iterations 
without the use of constant time steps fail to converge as well. 

3.2 Sharp, Double-Cone 

This unexpected result in which the ability to achieve a steady state requires 
temporal accuracy has not been observed previously in typical, convex blunt body 
geometries. However, in a recent hypersonic code validation study 1 , simulated flow 
over a shaip, double cone (25° / 55°) at Mach 9.5 and Reynolds number 264000 m" 1 in 
nitrogen showed large scale, periodic shedding of vortices in the separation region 
around the 30° compression. A 1 st order accurate (in space), constant Courant number, 
simulation showed a steady flow for this same case. A 2 nd order simulation at a lower 
Reynolds number (144000 m' 1 ) showed a steady flow result on the same grid. 
Experimental data" in the Calspan - University at Buffalo Research Center (CUBRC) 
Large Energy National Shock (LENS) tunnel for this case (Run 24) from thin-film heat 
transfer gauges indicated that the flow achieved steady state in the test time. The non- 
linear minmod flux limiter 3 in LAURA was a suspected source of numerical ringing 
within the recirculation region that triggered the unsteadiness. 





Ms 


Ml 






Figure 2: Pressure and heating over sharp, double cone for Run 24 of Reference 2 using time-accurate 
simulation to attain steady state. 

In light of these recent results showing that attainment of a steady state can depend 
on temporal accuracy in some circumstances, the shaip, double- cone case was revisited 
using time-accurate simulation. The time step was set at 3.65 10" 7 s. The reference flow 
time, L / V„ , is 7.08 1 (f s. The flow converged to a nearly steady state - very slight 
oscillation of streamlines persisted in the recirculation region, but these oscillations 
were orders of magnitude less severe than observed in the constant Courant number 
simulation. Comparisons with experimental data for pressure and heating are presented 
in Figure 2. Good agreement with the onset of separation is observed. Fair agreement 
with reflected wave peak pressure and heating location is obtained. The broad, flat 
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pressure plateau in the recirculation region for 0.5<x/L<l had significant, varying 
structure in the previous, constant Courant number simulation. 

Both the chute and double-cone cases have a subsonic region contained in a 
concavity (rays perpendicular to the surface converge in a subsonic domain) driven by 
an external supersonic flow. In the typical grid that is highly stretched from the body a 
constant Courant number simulation means that physical time advances more rapidly 
away from the body than in close to the body. A similar argument is suggested for cases 
in which a constant time step is used but temporal errors are largest in the near wall 
region if an insufficient number of sub-iterations are applied. Perturbations associated 
with each relaxation step may tend to aggregate near the wall in the concavity before 
exiting a supersonic outflow boundary. Waves accumulate (in this non-physical 
temporal map ), strengthen, and reflect back in to the subsonic domain in a self- 
sustaining but non-physical manner. This scenario is consistent with observed behavior 
in the two systems tested here. Certainly, not all subsonic domains with concavities 
induce this pseudo-unsteady behavior in constant Courant number simulations - there 
are many examples of compression comers with separated flow being computed with 
constant Courant number relaxation. It is simply noted that such flow conditions may be 
susceptible to pseudo-unsteady behavior in non-time-accurate relaxation and the 
overhead associated with implicit, time-accurate simulation may be recovered with 
better convergence. 

4 CONCLUDING REMARKS 

Point- and line-implicit relaxation algorithms within Program LAURA are reviewed. 
Modifications involving sub-iterations to achieve temporal accuracy are defined for the 
purpose of investigating their effect on convergence. Hypersonic flows with subsonic 
regions bounded by a local concavity are observed to support periodic, unsteady 
behavior using a constant Courant number relaxation but recover a steady flow behavior 
using time-accurate simulation. The test problems require a time-accurate simulation to 
evolve to a steady state. Relaxation using a constant time step but insufficient sub- 
iterations to attain first-order time accuracy fails to produce a steady state. Multiple sub- 
iterations without the use of constant time steps fail to converge as well. The behavior is 
thought to emerge from the temporal mapping of local time stepping allowing 
accumulated perturbations to reflect back into the subsonic domain. 
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